% Calculation for covariance matrix
function [C,Z] = CovMatriC(params)
% j_max: maximum terms of Zernike polynomials
% R = Radius for phase screen
% L0 : Outer scale [m]
% l0 : Inner scale [m]

j_max = params.j_max;
r0 = params.r0;
R = params.D/2;
L0 = params.L0;
l0 = params.l0;

k0 = 1/L0;
% indexs saving matrix generation
% ----------------------------
Z=zeros(3,j_max);                   %
for j=1:j_max
    idx = j;
    Z(1,idx)=j;                       % First row for j
    [n,m] = Noll_j_to_nm(j);
    Z(2,idx) = n; % Second row for n
    Z(3,idx) = m; % Third row for m
end

%------computing covariancce matrix C--------
%--------------------------------------------
C = zeros(j_max-1,j_max-1);
if (L0==inf) && (l0==0)
    % kolmorogrov spectrum
    for q=1:(j_max-1)
        for p=1:(j_max-1)
            if Z(3,p+1)==Z(3,q+1) && (rem(abs(p-q),2)==0 || (Z(3,p+1)+Z(3,q+1)==0))
                num=2.2465*(-1)^(Z(2,p+1)/2+Z(2,q+1)/2-Z(3,p+1))*sqrt((Z(2,p+1)+1)*(Z(2,q+1)+1))...
                    *gamma(Z(2,p+1)/2+Z(2,q+1)/2-5/6);
                den=gamma(Z(2,p+1)/2-Z(2,q+1)/2+17/6)*gamma(Z(2,q+1)/2-Z(2,p+1)/2+17/6)...
                    *gamma(Z(2,p+1)/2+Z(2,q+1)/2+23/6);
                C(p,q)=num/den;
            end
        end
    end
    C = C*(2*R/r0)^(5/3);
    
else
    % von-karman spectrum
    for j1 = 2:j_max
        idx1 = j1-1;
        for j2 = 2:j1
            idx2 = j2-1;
            if (Z(3,j1)==Z(3,j2)) && (rem(abs(j1-j2),2)==0 || (Z(3,j1)+Z(3,j2)==0))
                n1 = Z(2,j1); n2 = Z(2,j2); m = Z(3,j1);
                fun = @(k) (besselj(n1+1,2*pi*k).*besselj(n2+1,2*pi*k)./k)...
                    .*exp(-(2*pi*k*l0/(5.92*R)).^2)./(k.^2+(R*k0)^2).^(11/6);
                C(idx1,idx2) = (0.0458/pi)*(R/r0)^(5/3)*sqrt((n1+1)*(n2+1))...
                    *(-1)^((n1+n2-2*m)/2)*integral(fun,0,inf, 'ArrayValued', true);
                if j1 ~= j2
                    C(idx2,idx1) = C(idx1,idx2);
                end
            end
        end
    end 
end

% end of function
end


